arXiv:1505.03785vl [q-bio.CB] 14 May 2015 


A dynamical model of the adaptive immune system: 
effects of cells promiscuity, antigens and B-B interactions 


Silvia Bartolucci^ and Alessia Annibale 

^ Department of Mathematics, King’s College London, The Strand, London WC2R 2LS, UK 
^ Institute for Mathematical and Molecular Biomedicine, King’s College London, Hodgkin 
Building, London SEl lUL, UK 

Abstract. We analyse a minimal model for the primary response in the adaptive immune system 
comprising three different players: antigens, T and B cells. We assume B-T interactions to 
be diluted and sampled locally from heterogeneous degree distributions, which mimic B cells 
receptors’ promiscuity. We derive dynamical equations for the order parameters quantifying the B 
cells activation and study the nature and stability of the stationary solutions using linear stability 
analysis and Monte Carlo simulations. The system’s behaviour is studied in different scaling 
regimes of the number of B cells, dilution in the interactions and number of antigens. Our analysis 
shows that: (i) B cells activation depends on the number of receptors in such a way that cells 
with an insufficient number of triggered receptors cannot be activated; (ii) idiotypic (i.e. B-B) 
interactions enhance parallel activation of multiple clones, improving the system’s ability to fight 
different pathogens in parallel; (iii) the higher the fraction of antigens within the host the harder 
is for the system to sustain parallel signalling to B cells, crucial for the homeostatic control of cell 
numbers. 
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1. Introduction 


The immune system is a complex collection of organs, tissues and cells which is present in all 
vertebrates and protects the organism from external pathogens [T]. In this work we introduce 
a minimal model to describe the primary response of the adaptive immune system, a network 
of highly specialised cells that produces a targeted reaction to specihc antigens, i.e. viruses. 
The main players are B and T cells, two different types of lymphocytes. They independently 
recognise the antigen binding it with their receptors (£g. [^. Each group of T and B cells sharing 
the same receptors’ structure (clone) is able to recognise and hght only a particular virus with 
complementary epitope. 

The immune response by B cells is activated or suppressed according to a conhrmation signal 
sent by T cells in the form of excitatory or inhibitory proteins, the cytokines. Such response 
(when activated) consists in the secretion of antibodies, proteins able to chemically bind and 
neutralise the antigen, hence (possibly) avoiding the propagation of the infection. This two- 
signals mechanism prevents erroneous B cells activations against mismatching antigens or other 
cells of the organism. 

In recent years, this system composed by a large number of interacting agents - T, B cells 
and antigens - has been looked at through the prism of statistical mechanics to understand its 
global features and functionalities [a El HE]. Following this promising line of research we propose 
here a model for the dynamics of T and B cells and antigens, extending preliminary proposals (see 
m and references therein) to incorporate important biological features of real immune systems. 
In particular, we study the effect of having B cells with a variable number of receptors on their 
surfaces, one of the most important mechanisms preventing autoimmune reactions and diseases: 
a receptor editing process is indeed very commonly observed during the B cell maturation, where 
self-reactive cells, which may be responsible for the onset of autoimmune responses, are suppressed 
at an early stage of the development by changing the number of receptors [7|. We also introduce 
interactions between B cells, the so-called idiotypic network, hrst hypothesised by Jerne [8]: B cells 
receptors can not only recognise antigens but also other B cells with complementary receptors. 
In the presence of an antigen antibodies with complementary receptors, called idiotypes and 
denoted with Abl, are produced and recognised by complementary antibodies (Ab2) which share 
structural features with the original antigen. According to recent experimental studies PEnmu 
the idiotypic network seems to play a central role in autoimmune diseases, supporting a cascade 
of autoantibodies production, which recognise each other and modulate the immune response. 
Finally, we incorporate the effect of an external antigenic held to investigate the immunological 
memory [12] , i. e. the ability of the immune system to produce a more effective and faster response 
at a second encounter with an antigen. 

From a statistical mechanics perspective this system is modelled as a bipartite network, where 
links between B and T-cells are sparse as biologically required. We introduce a phenomenological 
Hamiltonian description of the system and carry out a dynamical analysis of the network, evolving 
via a Glauber sequential update. Via non equilibrium statistical mechanical techniques initially 
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Figure 1. Schematic representation of the antigen recognition process and immune response 
activation by T and B cells. The best-matching B and T cells independently detect the antigen: 
the active T clone sends excitatory cytokines (green links) to the B clone to start the antibodies 
production and suppress the non-complementary B clones with inhibitory signals (red). The other 
T cells are in a quiescent state. 


developed for Hopfield neural networks and spin glasses [ISl[Il[I51[ISl[n] we derive the equations 
for the time evolution of a set of parameters, quantifying the immune response strength. We 
analyse the system’s behaviour in different regions of the parameter’s space (he. number of 
clones and triggered receptors, noise level, etc.) through linear stability analysis and Monte Carlo 
simulations. The paper is structured as follows: in Sec. |^we dehne the model, in Sec. |^we give 
an overview of the results that we derive in Sec. H |5]and§ We summarise our conclusions inj^ 


2. The model 


We model T-cells as binary variables or “spins” cTj = ±1, i = l,...,iV, which can be active 
he. secreting cytokines (+1) or quiescent (—1). B cells and antigens are characterised via their 
concentration, respectively 6^ G K, /i = 1,..., P and 'i/’a G K, a = 1,... ,A with respect to a 
reference level. Both the number of B clones and the number of antigens are sub-linear in the 
system size N, with P = N^, 6 G [0,1) and A = iV“, a G [0,1). In the absence of antigens and 
interactions with T cells, the B clone sizes 6^ can be regarded as Gaussian variables na. Without 
loss of generality we take their mean to be zero and we denote by their covariance matrix. 

The interactions between the Pth T clone and the /r-th B clone, mediated via the cytokines, are 
represented by the variables G {+1, —1, 0}, respectively corresponding to excitatory, inhibitory 
or absent signals. We will regard those as random variables drawn from 


^(«)=n 

2/i 


% 

2N'y 





%,o 


= 0(iV°) V/r, (1) 


where the g^’s are drawn from a distribution V{q) = P control the degree of B 

cells promiscuity, he. their ability to communicate with different T cells, via different receptors. 
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The fraction of non-zero B-T links determines the degree of dilution of the system: for 7 = 0 
the system is hnitely diluted, whereas for 7 > 0 the system is extremely diluted. The combined 
interacting system of B-T cells and the antigens can be phenomenologically described by the 
following Hamiltonian: 

p / A N \ p 

b) = - ^ 6 J ^ ^ ]+—J2 b^A^K. (2) 

fj,=l Va=l i=l ) ^ 

The hrst term takes into account the interactions between antigens and B cells via the matrix 
7 ^, the second term is related to B - T interactions mediated by cytokines, the third one takes 
into account the effect of the idiotypic network (B-B interactions) (£g. [^. B-B interactions are 
mediated via A = {A^y}\ according to the theory of idiotypic interactions [HI EO], B clones can 
recognise not only antigens but also antibodies with complementary epitopes creating a network 
of imitative interactions between B cells. We represent epitopes as binary strings and assume 
that complementary strings, like e.g. 010 ... and 101 ..., excite each other. Also, we assume that 
we can order the strings on a ring in such a way that each string sits close to affine strings and 
opposite to complementary ones (£g. right). Hence, we suppose that the /x-th B clone expansion 
is triggered by the (/i -|- P/2)-th B clone, which is precisely complementary to that B clone (£g. 

left). While complementary B cells excite each other, we assume that each B cell suppresses 
itself, to prevent uncontrolled production of a single cell type. Therefore, we are led to use for the 
B-B interactions the Toeplitz matrix 


k6, 


flU 


mod P? 


with k G [0,1) representing the strength of idiotypic interactions. Its inverse is 

1 . k 

'C^,(i/+P/2) mod P- 


(A 


- 1 > 




+ 


(3) 


(4) 


' 1 -P 

We note that A is positive dehnite and symmetric and the following relations hold: A^y = A{fi — u) 
with A{n + P) = A{n) and A{n) = A{—n). The eigenvalue spectrum p{a) of A has a hnite limit 
as N ^ 00 , ensuring the correct scaling of the Hamiltonian. B-antigen interactions are taken 
into account via the matrix 7 ^: antigens will excite complementary B and they will repress the 
identical one, in a similar way as B clones detect and excite each other as shown in £g. We 
can further investigate the B-Antigen interactions analysing the stochastic process that governs 
the dynamics of B cells concentration. Assuming that the B cells concentration evolves according 
to a gradient descent on the Hamiltonian ([^, we have 

A N , P 


dt 


■ — + noise = ^ 

^ a=l i=l ^ ^ v=l 


(5) 


Denoting with 7 the concentration of the antigen complementary to the clone, we need the 
concentration to increase when 7 ^ 0 . Also, we assume that clone 7 +p/ 2 , complementary to 7 
and thus carrying the same epitope as antigen is inhibited by the presence of antigen 7) 
we denote ki the strength of this inhibition. This leads to a matrix for the B-Antigen interaction 
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in the form 


'IJa — ^tJ.a ^l^/i,(a+P/2) mod P 


( 6 ) 


According to (|^, the concentration of the /x-th clone also increases in the presence of excitatory 
signals received by T cells (second term) while the third suppressive term is related to B-B 
interactions acting as a threshold to be overcome to start the immune response. Assuming that 
the total B clones concentration is conserved on average leads to a relation between k and ki 


dt 




= 0 ^ 


1-ki ^ 1 

l-k 2y/P '>Pa 


(7) 


which depends on the steady state concentrations of B cells and antigens and on the inverse noise 
level For simplicity we will set k = ki, which leads to r/ = A and to an equilibrium ratio 
between B cells and antigens concentration only controlled by the noise level in the system. 


B cells T cells 


Antigen 





Figure 2. Schematic interactions between B, T cells and the antigen A. In the presence of an 
antigen with concentration ijj, the complementary B cell will detect it (B-A interactions mediated 
via the matrix t]) and will receive a confirmatory signal from the active T cells (represented by up 
arrows) via the cytokines . 


2.1. Dynamical equations 


At equilibrium at inverse noise level (consistent with our assumption b ~ 7V(0, A ^)) we 
expect the joint distribution P{cr, b) to be given by the Boltzmann distribution 

-^H[cr,h) 

Pi^,h) = -^- ( 8 ) 

The equilibrium marginal distribution for the cr is found, by integrating out the variables 6^, 




1 /3 


E,, E^,. Ei. E^,. 


( 9 ) 
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Figure 3. Left: B-B interaction via Af^^. The expansion of the /i-th clone is triggered by its 
complementary clone. Right: We represent the epitopes as binary strings, organising them on a 
ring and assuming that the complementary strings interact. 
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Figure 4. Scheme of B-B and B-A interactions. Antigens and B cells are denoted by variables 
0,1 representing the shape of their receptors: 0—1 variables represent complementary receptors 
(key-lock mechanism). Different B cells excite each other and each of them represses itself, while 
the antigen will excite complementary B repressing the identical one (0-0). 


in the Boltzmann form with effective Hamiltonian 'H{cr), involving only interactions between T 
cells, 

W(«j) = X) Y, + A'’-' E (10) 

ij ia 

with the separable form Jij = where and thus describing 

an associative network with diluted patterns encoding hghting strategies against different 

antigens |6]. In the regime we consider here 5 < 1, the associative network is away from saturation. 
Analysis near saturation have been carried out in statics for A = 1, mostly for V{q) = S{q — c) 
IIH1II9]. We can rewrite the Hamiltonian ([l0|) as 


'H(cr) = —-M^(cr)A ^M(cr)-|-^M(cr) 


( 11 ) 
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(12) 


in terms of the order parameters M = (Mi,..., Mp), where 

1 ^ 
i=l 

quantihes the strength of the excitatory signal on B clone /i and thus its activation. Here 
t/? = {'ipi,... ,'ipp) and we denotes the transpose of v. Next, we assume a Glauber 

sequential dynamics for the variables cr, converging to the equilibrium measure P{cr), so that 
the instantaneous probability of hnding the system in state cr = (ai,..., a^) at time t is governed 
by the master equation 

- [Pt{Fi(T)Wi{Fi(T) - Pt{(T)Wi{(T)] , (13) 

2=1 

where F* is the Gth spin-flip operator Fj((Ji,Uj,..., a^) = (ci, ■■■, —o'i, ■■■cr^) and transition rates 
between cr and Ficr have the Glauber form 

wM = ^[1 - aitanh{/3hf{(T))] , (14) 

with effective held 


Af = iV’-i (15) 

ensuring convergence to (§. Here T = 1/13 represents the rate of spontaneous spin-hip, hence the 
ehective noise of the T cell dynamics. Using ( |I^ , we write the dynamics for the B cells activation, 
described by the probability of hnding the system in a macroscopic state M at time f, namely 


P,(M) = ^ Pi(o-)5(M-M(o-)). (16) 

cr 

Via a Kramers-Moyal expansion for large system size and away from the saturation regime it is 
possible to show that Pt(M) evolves according to a Liouville equation [ 6 ]. It follows that M evolve 
deterministically with a dynamics described by 


dM 

dt 


= ( ^^tanh [(3 ) - M 


(17) 


where C = r/A ^ and ( ■ )^ denotes the average over the distribution (1). For rj dehned in ([^ 
and k = ki, we have C = 1, and each order parameter evolves accordin 


dM^ 

dt 


= ( tanh 


(3{M^A-^^ + 


M„ 


(18) 


Diherent choices of the matrix 77 or of the constants k, ki will lead to diherent forms of the matrix 
C, however the choice C = 1 seems to ensure the strongest excitatory signal to antigen-activated 
B clones. As an illustration, let us consider the simple case where we have just one antigen i.e. 
'ipi'. for the choice C = 1 the held acting on the i-th T cell exciting the complementary clone hi is 
'ipi, hence i receives an activation signal only from the clone B activated by the antigen. With two 
antigens present, whose complementary B cells are signalled by diherent cytokines from the same 
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Figure 5. Left: effect of having C = 1, in the presence of two antigens whose complementary 
B cells are signalled by different cytokines from the same T cell: the field acting on the T cell 
is tpi — 'ip 2 and T gets activated only if ipi > 'ip 2 - Right: an illustration of the effect of having 
C = which may result in an inhibitory signal on T cells from inactive B cells. 


T cell, the field acting on the latter is -01 — '02 and the T cell will get activated only if ipi > 02, as 
biologically desired (fig. left). For C 0 1 there might be a negative interference on the desired 
signal. For example, for the alternative choice t] = 1, leading to C = A“^, the field acting on T 
in the presence of a single virns 0i 0 0 wonld be meaning that for = +1, 

= — 1 the T cell would receive, besides the excitatory signal from the B cell activated by the 
antigen, an inhibitory signal from a non-activated B-cell (fig. right). As a result the overall 
field and immune response will decrease. Similar arguments can be given for rj defined in (|^ with 
ki 0 k. 


3. Overview of results 

The sparsity of the B-T connections makes the system able to activate multiple B clones in 
parallel. This multitasking capability is one of the core features of the immune system that, 
in normal conditions, can control and block several simultaneous antigenic invasions. We find 
that the parallel activation of B clones may occur in symmetric fashion, where all infections are 
fought with the same strength, or in a hierarchical fashion, where the system prioritises immune 
responses against specific pathogens. We are able to identify the system’s parameters, such as 
noise, number of links, number of different infections, which induce the switch from the symmetric 
to the hierarchical operational mode. This might be potentially useful to investigate the causes of 
dramatic failures in the functioning of the immune system. The switch to a hierarchical immune 
response is in fact mostly related to the presence of strong infections, such as hepatitis and 
autoimmune disorders: while the immune system invests the highest amount of resources tackling 
the main disease, the progression of minor infections may become lethal [2T1122 ] . 

An important feature of the model is the dependence of the B cells activation on the number of 
receptors on their surface. Results in Sec. |^show that cells with few receptors may be transiently 
activated but fail to sustain the signal even if strongly excitatory (fig. [^. The number of triggered 
receptors affects both the immune response strength and the critical temperature at which they 
get activated: both decrease with the number of triggered receptors. In addition, competition 
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emerges between B clones to be activated. As the number of activated B clones increases, signalling 
pathways to inactive clones get more noisy due to the interference of active clones. The critical 
temperature for the activation of inactive clones is a decreasing function of the fraction of active 


clones (73). 


Another result of biological interest is that idiotypic i.e. B-B interactions contribute to the 
overall stability of the immune system, preventing unwanted activation and increasing the region 
where all clones are equally activated and ready to start an immune response upon arrival of 
new infections (Sec. [^. In particular, including B-B interactions in the model affects the critical 
temperature, in this case widening the region where a symmetric immune response is stable, as 


the interactions strength k is increased (£g. 22). 

Finally, the role of antigens in B cells activation is understood by our model as that of 
external helds on coupled ferromagnetic systems One of the interesting consequences is then 
the presence of hysteresis phenomena [23] , which could explain short-term memory effect in the 
immune response HH EH, even in the absence of memory cells. The effect of antigens on the 
immune system basal activity and surveillance is also of interest: the response of non-infected B 
clones decreases with the fraction of infected clones, due to the interference of strongly activated 
B cells, and is activated at a lower T, making the whole system unresponsive to new incoming 


viruses (£g. 23). 


4. Effects of receptors’ promiscuity 

We hrst consider the case where there is no antigen = 0 V a, and no B-B interactions i.e. A = 1 
and focus on the effect of having a variable number of triggered receptors on different B clones, 
i.e. heterogeneous in ([^. In order to compare activation of B clones with different numbers 
of receptors, it is convenient to look at the activation per receptor, given by the normalised order 
parameters /i = 1,... ,P, which take values in the range [—1,1] for all clones fi. 


The dynamical equations (18) then read 


dm^ 

dt 


N'T / 

= — (^^/.tanh I 



m, 


(19) 


One can show that at the critical temperature Tc = gmax, where gmax = niax^[g^], the system 
undergoes a phase transition, with the equilibrium phase at T > Tc characterised by m = 0, and 
B clones activation occurring at T < Tc, where m 7^ 0. At the steady state (dm/df = 0), we have 

M, = A^(etanh(/?5^rM.))^ . (20) 

U 

Taking the scalar product with M and using the inequality | tanha;| < |a:|, yields 




< 




( 21 ) 

U fl 

This implies = 0, hence m = 0, for /dgmax < 1, meaning that none of the B cells can get 
activated for noise levels above the critical value T > gmax- Although m = 0 is a steady state 
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solution of (19) for any value of T, a linear stability analysis shows that it becomes unstable for 
T < gmax- To this purpose, we compute the Jacobian of the linearised dynamics about the steady 
state m* = 0 


(m) 


which gives 


J jjLU 


drrii 

% 


N"/ 

= — (^'^tanh(/3^^'^g‘^m^))^ - (22) 

^ U=1 




(23) 


Substituting m* = 0 we get J^u = {(3q^ — l)6fj,y. This is a diagonal matrix, hence the largest 
eigenvalue, which gives the stability of m* = 0, is Amax = /^Q'max — 1- This gets positive for 
Pqmax > 1, showing that non-zero solutions m 7 ^ 0 will bifurcate away from m = 0 at T = gmax- 
We inspect the structure and the stability of the bifurcating solutions hrst for a toy model with 
just two B-clones and then for the general case with P B-clones. 


4-1. A toy model with two B-clones 


Here we study a toy model with P = 2 B clones, and assume qi> q 2 - For 7 > 0 this model reduces 
to two independent Curie-Weiss ferromagnets, with critical temperatures qi and q 2 respectively 
(see discussion in Sec. 4.2). Hence, the most interesting case is obtained for 7 = 0 . The state 
m* = ( 0 , 0 ) is the only steady state of the system for T > gi, but it destabilises for T < gi. 
We can understand the system’s behaviour below criticality by solving numerically its dynamical 
equations 


dmi 

dt 

dm2 

dt 


0.2 

(1 - g 2 ) tanh(/9gimi) — [tanh(/3(gimi -h q 2 m 2 )) + tanh(/?(gimi - g 2 m 2 ))] - mi , (24) 
(1 - gi) tanh(/dm 2 ) -h Y[tanh(/3(gimi g 2 m 2 )) - tanh(/?(gimi - g 2 m 2 ))] - m 2 . (25) 


In £g. 1^ we show the flow diagram and the stable hxed points at different temperatures: first we 
notice that B cells with a higher promiscuity produce a higher immune response (mi), whereas 
a lower promiscuity g^ results in a lower or null activation (m 2 ) depending on the temperature. 
Hence, the number of receptors on B cells’ surface affects the responsiveness of B clones. In 
particular, if T is high, only clones with the highest number of receptors are activated and 
the system’s hxed point corresponds to the pure state mi 7 ^ 0, m 2 = 0. Lowering T induces 
the activation of cells with fewer receptors but with a lower intensity (mi > m 2 ). Theoretical 
results are consistent with Monte Carlo simulations, shown in £g. Furthermore, Monte Carlo 
simulations match experimental results [25] showing that cells with very few receptors are triggered 
transiently but fail to be activated in the long run, since the number of receptors is not sufficient 
to allow these cells to sustain the signal, even if strongly excitatory. Fig. shows that our model 
can reproduce this effect: cells with few receptors (green) have a lower activation even if triggered 
by a strong signal (initial condition) and tends to be switched off after a short transient, conversely 
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cells with a higher number of receptors (blue) produce a strong immune response, even if triggered 
by a low signal (initial condition). Analytically, we can investigate the structure of the hrst states 



- 1.0 - 0.5 0.0 0.5 1.0 

m, q, 



- 1.0 - 0.5 0.0 0.5 1.0 

m, (j. 


Figure 6. Phase portrait of the dynamical system (p4{ 25) for qi = 0.6, q 2 = 0.4 at low 


temperature T = 0.01 (left) and high temperature T = 0.2 (right). Red lines represent null-clines 
and stationary states are at the intersection of null-clines. 




Figure 7. Monte Carlo simulations with 10^ spins for qi = 0.6, <72 = 0.4 at low temperature 
T = 0.01 (left) and at high temperature T = 0.2 (right). Markers represent the numerical solutions 

of @1^ . 


to bifurcate below Tc by expanding the steady state equations, obtained setting dmi^ 2 /dt = 0 in 
(24, 25), for small mi, m 2 , close to criticality, i.e. at /3gi = 1 -|- e. We get as possible solutions 


m\ = 3e — 3^7712 

ql 


rUn 




= 3^(l + e) 

H2 


4 

^2 


2 

3^mi 

T2 


(26) 


(27) 


as well as mi = 0 ,m 2 = 0. Solution (27) is unphysical as it stays 0{1) at e = 0 for qi 7 ^ q 2 , hence 
m 2 = 0. Inserting m 2 = 0 in (26), we get ml = 3e. Hence, the hrst state bifurcating away from 
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Figure 8. Monte Carlo simulations with N = 10^ spins, with qi = 0.6, q 2 = 0.4 at T = 0.2. Clone 
activation mi, m 2 as a function of time for different initial conditions. Cells with few receptors, 
q 2 (green) fail to be activated even if triggered with a strong signal. 


m = (0,0) is in the form m = (\/3e, 0). At high temperature (below criticality) only the clone 
with the higher number of receptors is switched on. For qi = q 2 we clearly retrieve the results 
in [6] and both B clones are activated with the same intensity below criticality. Next we derive 
the region in the phase diagram where cells with fewer receptors become responsive. Naively one 
might expect that m 2 becomes active at T ~ g 2 - In reality, heterogenities in the cell promiscuities 
deeply affects cell responsiveness and cells with fewer receptors will remain quiescent even at very 
low T. By Taylor expanding (19) for small m 2 in powers of e = q2{3 — 1, we obtain 

I (28) 

(29) 


mi = tanh(/dgimi) ~ tanh ( —mi 

\<12 


m: 


2 = 3e - 


Since mi is 0{1) the non-zero solution for m 2 is impossible close to T = q 2 . The pure state 
m = (mi, 0) will then have a wider stability region, that can be found by analysing the eigenvalues 


of the Jacobian (22) at m* = (mi,0). 


-^1 = PQi ~ f3qitanh^{f3qimi) — 1 , 

(30) 

A 2 = l3q2 - f3qiq2 tanh^{j3qimi) - 1 . 

(31) 

Analytically we can calculate Ai ^2 near T ~ 0 


Ai ~ -1 , 

(32) 

A 2 ^ (3q2 - (3qiq2 - 1 • 

(33) 

and near Tc, i. e. at Pqi = 1 -|- e 


Ai ~ -2e , 

(34) 

^ ® - 91 

(35) 


Qi 
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For intermediate T we compute Ai ^2 numerically. Plots of eigenvalues as a function of the 



Figure 9. Eigenvalues (31) as a function of T for qi = 0.8 and q 2 = 0.7. Left: Ai, the red dashed 
line represents the behaviour near T ^Tc (35). Right: A 2 T, the red marker represents the limit 
at T ~ 0 dMl). 


temperature are shown in hg. where the theoretical predictions for T ~ 0 (33) and T ~ Tc 
(35) are highlighted. We note that Ai < 0, V T, hence the stability of the pure state is determined 
by the sign of A 2 . In £g. [T^we show a contour plot of A 2 = 0 in the T — q 2 plane for different 


values of qi. The linear behaviour can be understood as follows. In the pure state region we have, 



Figure 10. Phase diagram in the space (T, 92 ) fixing qi = 0.8 obtained from the condition A 2 < 0 
(31). The dotted line represents the theoretical critical temperature (38). In the (PS) region the 
pure state is stable (mi ^ 0 ,m 2 = 0). At low T, B clones are hierarchically activated (H). 


using the steady state equation mi = tanh(/3gimi), 

Ai = /?gi(l - mi) - 1 , (36) 

A 2 = /^g2(l - gi^i) - 1 . (37) 

As T decreases (below ^ 2 ) 'mi increases, so that the eigenvalues stay negative and stability is 
ensured, until mi reaches its maximum value mi = 1. At this point, a further decrease of the 
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Figure 11. Plot of mi,m 2 as a function of < 72 . Left : qi = 0.8 at T = 52 , Right: <71 = 0.8, 
T = 0.01,0.05,0.1,0.3. 


temperature will make A 2 positive, destabilising the pure state, and a new state (mi, m 2 ) will take 
over. The temperature at which bifurcations from the pure state are expected is thus found from 
the condition A 2 (T, mi = 1) = 0, as 


T = ®(1 -?i) 


(38) 


which is in agreement with the critical temperature computed numerically (£g. 10). Deviations 
from the linear behaviour are expected in the regime qi — q 2 , where a symmetric activation occurs 
for gi < 1/3 [6]. In £g. 11 we plot mi, m 2 as a function of q 2 in the (PS) region (left) and crossing 
the critical line where m 2 becomes non-zero (right). In the latter region, i.e. for q 2 > T/{1 — qi), 
the stable state is m = (1, m 2 ) where m 2 is the T-dependent stationary solution of ( |25| ) at mi = 1. 
In particular, for T = 0 the stable state is 


mi = 1 , (39) 

m 2 = 1 - , (40) 

in agreement with simulations and flow diagrams (fig. |^[^. In conclusion, the system can activate 
clones with different numbers of receptors simultaneously for T < ^ 2(1 — qi)- The activation is 
hierarchical (H), with clones with higher promiscuity being prioritised with respect to the others. 
In particular, clones with the highest number of receptors are activated with the strongest possible 
signal in a wide region of the phase diagram. 


4-2. The case of P B clones with a variable number of receptors 

Next we study the case where the number of B clones is P = where 6 G [0,1), and N is the 
number of T-clones. The dynamical equations are 

^ tanh - m^ , (41) 

which can be rewritten as 

tanh [ I3{q^m^ + Y ^ • (42) 
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Upon introducing the noise distribution Pf^{z\{m^,q^}) = {6{z — ^ clone p, we 

have 

QT?7 /* 

= I dzPi,{z\{m^,}) tanh + z)) - (43) 

where Pfj,{z\{m^,q^}) can be written, using the Fourier representation of the Dirac delta and 
carrying out the average over as 


/ + 00 J 


/ + 0O J ^ 

.(44) 

-oo 27r 

If 7 > 0, extending the sum in the exponent to all patterns will add a negligible contribution in 
the thermodynamic limit, hence, as iV —)■ oo all clones will have the same noise distribution 

p 


P^l{.A{mu,qy}) Pp(^|m,q) = {6{z -^q^Cm^))^ 


(45) 


U=1 


We note that the sum on the RHS of (45) is at most Cl(~ '^), hence for 5 < 7 it is negligible 


in the thermodynamic limit. This yields Pp( 2 ;|m, q) —)■ 6{z) as iV —)■ 00 , so that equations (43) 


decouple and the system reduces, for 5 < 7 , to a set of indipendent Curie-Weiss ferromagnets, 
each evolving according 
cIttT/ 

= tanh(/3g^m^) - m^. (46) 

At the steady state, each B clone fi becomes active at its own critical temperature = q^, 


independently of the other clones (fig. 12, left). In contrast, for 5 > 7 the noise distribution 
Fp(z|m, q) has a finite width, due to clonal interference, and the equations for the evolution of 
clonal activations are coupled. In this regime, clones compete to be activated and the ones with 


fewer triggered receptors will fail to be switched on (£g. 12, right). In the following section we 


will analyse the effect of receptors heterogeneity in the regime of competing clones. We will show 
that clonal interference affects both critical temperature and intensity of B clones activations. 


J 1 .. 2 .I. Bifurcations near the critical temperature and stability region in the regime of competing 
clones. In this section we study the bifurcations away from m = (0,..., 0) below the critical 
temperature Tc = gmax- Without loss of generality, we assume gmax = gi- We Taylor expand the 
steady state equations obtained by setting dm^/dt = 0 in (41), for small at /3qi = 1 -|- e 






m„ 


% 


(3^ 

% 3 




upX 


0 P^qfi 3 n qfi 2 

= Pq^rup - —m^ - ^ 


2 3 


For any /i, = 0 is always a solution. Non-zero solutions are given, for /i = 1, by 




(3qi-l- 


2 I l^qf 2 p 2 3 n 

p>i 


N-i' 


(47) 

(48) 

(49) 
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Figure 12. Monte Carlo simulations with N = 10'^ spins, q^j, = e^q. T = 0.01, e = 0.7, qi = 6. 
Left: 6 = 0 . 2,7 = 0.25. In the non-competing regime (i5 < 7 ) all clones are activated. Right: 
5 = 7 = 0.25. In this case only few clones with the highest number of receptors are active. 


that at /Sgi = 1 -h e and for large N yield 

P>1 

For yU 7 ^ 1, we have 

- 


2 

iV7 


P 2 3 n 

p>i 


(50) 


(51) 


which gives, for —)■ cx) and to C>(e°), 


K = - ql) - -J^^rnlq^ + 0{N . 

V p>l 


(52) 


Snmming over /r > 1 we get 


+ =31 5^(gp-gi)] (53) 

p>i ^ p>i ' ^ p>\ ^ 

Since q^ < qi, V/i > 1, this equality can never be satished, showing that = 0,V/i > 1. 
Substituting this result into ( [^ , we hnd ml = 3e, hence the hrst state to bifurcate away 
from m = (0,..., 0) is m = (\/3e, 0,..., 0). Next we calculate the region where the pure state 
(mi, 0,..., 0) stays stable, by inspecting the sign of the eigenvalues of the Jacobian of the linearised 
equations of motion about the steady state. 

For a steady state with the general structure m = (mi, m 2 ,..., m„, 0,..., 0), where n is the 
fraction of activated clones, the Jacobian (23) has a block structure, where diagonal terms are. 
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for yU < n 


<4^ = “ (tanh^(/?(gy,my, + '^CQv'^u)))^) - 1 

U^ll 


and for /i > n 


J'/.M = /5g/.(l - (tanh^(/?^^^g^m^))e) - 1, 


while off-diagonal elements are, for fi,TJ < n 

PQu 


jiiu 


iV7 


(tanh^ (/3 C(Iv'^u))p - 


(54) 


(55) 


(56) 


and = 0 otherwise. For iV —)■ cx) the matrix becomes diagonal with eigenvalues given 

by (54) for /i < n and (55) for fi > n. In the pure state m* = (mi, 0,..., 0), where only one clone 
is activated, we have 

Ai = /3qi{l - (tanh^(/3gimi))^) - 1 , 

“ (tanh^(/?^^gimi))^) - 1 /i > 1 . 

Near the critical temperature Tc = qi, setting (3qi = 1 -|- e and using m\ ~ 3e gives 


-^1 = — 1 — I3^rn\q\ = —2e < 0 




A/. = - 1 - /^ = (1 + e)— - 1, 


/i > 1 


(57) 

(58) 

(59) 

(60) 


showing that for < qi the pure state is stable near criticality, as opposed to the case = gi, V/i 
studied in [6j, where all clones are activated with the same intensity below Tc- In the opposite 
limit T —)■ 0, we get from (57), (58) 

Ai ~ -1 , 


\-PqA 1 - I “ ^ ^ ■ 


(61) 

(62) 


showing that the pure state is unstable at low temperature. 

Indeed, as the temperature is lowered below gi, we expect clones with fewer receptors to 
get active. In particular, at T = 0 we expect all clones to be activated, in a hierarchical 
fashion, whereby the system sends the highest possible signal to the clone with maximum number 
of receptors, the clone with second highest number of receptors is signalled by the remaining 
N — spare T cells and so on. This leads to the following heuristic rule for noiseless clone 

activations 

A=i.....p. (63) 


= 


i=i 


in agreement with Monte Carlo simulations shown in fig. [T^ 
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Figure 13. Monte Carlo simulations with N = 10^ spins, = e^q. 
6 = j = 0.25. The markers represents the theoretical predictions at T 


T = 0.001, 

= 0 (p). 


e 


0.7, qi = 6, 


4-2.2. Sequential B-clones activation: critical temperature and interference effects. In this section 
we calculate the critical temperature at which clones /i > 2 with fewer receptors get activated. 
We focus on the the regime of clonal interference 5 > 7, as for 5 < 7 each clone /i gets active at 
its own critical temperature p = Without loss of generality we can set P = aN^ for 5 > 7, 
where a = 1 for <5 = 7 and a —>■ 00 for 5 > 7. In the following we will consider qi > q 2 > ■ ■ ■ > qp. 


Assuming that m bifurcates continuously from the pure state, we can Taylor expand (41) at the 
steady state for small nit,, with 7^ 1, while mi = 0{1). For /r 7^ 1 we have 


m,, 


N'^ / ( S ^ 

= -( ^^tanh f3{f^qimi p'^qu^ 

\ \ u=2 

= /^^tanh (f^quCrr 



u=2 




+ 0{N-'^) 


(64) 


p p 

p^v 


v=2 


^ 3 f^q^rnp ^ 3 2 

^ --—— 2^ 


(65) 




The solutions are m„ = 0, or 


1 Q 2 z-' , 


P % 3 2 

p¥=v 


( 66 ) 


Hence, for 1 — fqp > 0, m^ = 0 while for 1 — fq^ < 0 the /r-th clone may be activated. Hence, the 
hrst state to bifurcate away from the pure state is m = (mi, m2, 0,..., 0) at T = ^2- Its amplitude 
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at I3q2 = 1 — e is 


m\ ~ 3e 


Ml 

qlN^- 


-mi 


(67) 


i.e. ml = 3e+0{N~'^). As T is lowered below q 2 we expect that the clones will activate sequentially 
one after the other, each one at its own temperature. In particular, assuming mi,2 = we 

have for yU > 2 

p 


= 


tanh {l3^qyi 



+ C>(At-T') (63) 

i/>3 / ' $,,^'"¥=0 

and expanding for m,^>2 small at T < ^2 shows that m3 becomes non-zero at T = ^3-1- 
Generalizing for n N'^ activated clones mi,..., m^ = 0{1) we have for /i > n 

p 


rrif, = 


tanh (5 q^i 


u>n 





+ 0{nN-^) 


(69) 


giving as bifurcation temperature T„ = qn+0{nN~'^). As the number of activated clones increases, 
their cumulative effect on the activation temperature of the remaining clones increases and can no 
longer be neglected for n = 0{N'^). The activation temperature of pattern n -|- 1, when n ~ N'^ 
clones have been activated, can be worked out from 


rrin+i = (|tanh ^/3(g„+im„+i + guj'"^^ , 

that gives, upon insertion of J dz6{z — X)"=i = 1 

m„+i = / dzPn{z\m,ci)ta:nh{(3{qn+imn+i +z)) . 


(70) 


(71) 


with Pn{z\'m.^ q) dehned in (45). Taylor expanding for m,i+i small we have, to leading orders 


m„+i = / dzPn{z\m,ci) 


(1 - iaidc?{(dz))l3qn+imn+i + C>(m^ , J 


(72) 


where we have used P„(^|m, q) = P,i(—2;|m, q). A solution is m^i+i = 0 and a non-zero solution 
is possible for 

Mn+i \ — j (\zPn{z\rn,c^tajdc?{(dz) ^ ^ 

For n -C N'^ 1 Pn(2;|m, q) = S{z), and we retrieve Mn+i = 1 for the temperature at which m„+i 
becomes non-zero. For P„(z|m, q) having a small but hnite width we can use tanh(/3z) ~ /3z 


(3qn+i ^ 1 + 


<ll+i 


~ 1 -|- 


E n 3 

u=i 






(74) 


showing that Tn+i < Qn+i and deviations from g„+i depend on the promiscuity distribution of 
the activated clones. Equation (73) shows that as more clones are activated, these create an 
interference, encoded in P„(z|m, q), that decreases the activation temperature of the inactive 
ones. Furthermore, it suggests that the number n of clones that the system can activate {i.e. the 
number of 0{1) order parameters m^) at small but hnite temperature is 0{N'^). 
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4-2.3. Numerical examples. In this section we test (73) and look at the effect of receptors’ 
heterogeneity on the intensity of B clones activation, for three simple cases that can be treated 
analytically, with qi > q 2 and P = aN'^: 

(i) ■■■ ,< 12 ). only one clone has a higher nnmber of receptor; 

(ii) q = {qi,... ,qi,q 2 ,... ,q 2 ). half of the clones have promiscnity qi, and half have promiscnity 

q2] 

(hi) q = {qi,... ,qi,q 2 ): only one clone has a smaller number of receptors. 

Our goal is to analyse the increasing interference effect due to the activated clones with more 
receptors on the quiescent ones with less receptors. According to (73) active clones should play 
the role of interference terms that lower the critical temperature of the quiescent clones. 

(i) Clones with the same number of receptors will be activated with the same intensity [6] and 


the activation vector bifurcating away from the pure state will have the form (see Sec. 5.2.1 


for a rigorous derivation) m = (mi, m2 ,..., m2). The amplitudes mi and m2 are found from 

p 


mi = 


N'y 

qi 


tanh 13 {^^qimi + ^2 ^ C 


u=2 


= ( tanh I f3{qimi + q 2 m 2 ^ C 

' \ u=2 

= '^Pp-i{z\q)ta.iah(3{qimi + q2m2z) 



(75) 


and 


m 2 = ^Pp_2(2:|q)tanh/3[gi^^mi + g2"i2(l + z)] 

Z 

^ ^ Pp- 2 {z\(l) tanh P[q 2 ‘m 2 {l + z)] + 


(76) 


where 




(77) 


Using (44), we can write Pp{z\q) = e ®“/z(q'2a) where Iz{x) is the modified Bessel function 
of the first kind m, and Pp_i(z|q) ~ Pp- 2 {z\q) ~ Pp{z\q). The activation temperature of 
the clones with smaller promiscuity follows from Taylor expansion of ( |76| ) for small m2 giving 

m2 = '^Pp-2{z\q)[l3q2m2{l + z) - ]-{(3q2m2f{l + zf] + 0{N-^) 


= (3q2m2 - -{I3q2m2f{l + 3(z^))] + 0{N 
= (3q2m2 - ^{I3q2m2f{l + 3Q;g2)] + 0{N~^) 


(78) 
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where we used the parity of Pn{z\q) = P„(—z|q) V n and {z^) = — ®5'2- We 

retrieve (5q2 = 1 + 0{N~'^) for the activation teperature of m 2 , consistently with (73), for 
n = 1. The activation intensity at f3q2 = 1 + e follows from (78) as 

3e 


m 2 = 


1 + 3^20’ 


(79) 


as one expects for homogeneous clones in the absence of clone /x = 1 [6]. In £g. 14 we plot 
the amplitudes mi,m 2 resulting from (75), (76), as a function of the temperature (left) and 
those resulting from Monte Carlo simulations, as a function of time (right). The latters are 
seen to relax to the theoretically predicted steady state. A contour plot of m 2 = 0 in the 




Figure 14. Left: steady state solutions mi,m 2 of the dynamical system (75), (76), as function 
of the temperature T, for qi = 0.8, q 2 = 0.4. Right: Monte Carlo simulations with N = 10^ spins, 
5 = 7 = 0.25, T = 0.3, qi = 0.8, <72 = 0.4. The markers are the theoretically predicted steady state 
activations (fz^, ([76|). 


T — q 2 plane, computed numerically from (76), is shown in £g. 15 Deviations from the line 


T = q 2 are consistent with finite size effects N~'^. 

In conclusion, the presence of one clone with a higher number of receptors, does not affect, 
in the thermodynamic limit, the activation temperature nor the activation intensity of clones 
with fewer receptors. 

(ii) In this case we expect a transition from the state m = (mi,..., mi, 0,..., 0) to m = 
(mi,... ,mi,m 2 ,... ,m 2 ), with 


mi = E -Pp/2(D|q)Pp/2(^2|q) tanh(/ 3 (gimi(l + Zi) + ^2^22:2)) , 

^ 1,^2 

m 2 = ^ Pp/ 2 (^i|q)Pp/ 2 (^ 2 |q) tanh(/3(g2m2(l + Z 2 ) + gimi^i)) , (80) 

7:1,22 


The temperature at which the on-set of non-zero m2 occurs is found by Taylor expanding 
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Figure 15. Contours of constant m 2 with q = {qi,q 2 ,..., q^)- The contour indicating the on-set 
of a non-zero value of m 2 is very close to the blue dotted line T = q 2 theoretically predicted. 
Deviations are within finite size effects 0{N~^). 


(80) for small m 2 with mi = 0(1) 

m 2 = ^ Pp/2{zi\q)Pp/2{z2\q) [ta.nh{(3qimiZi) + (1 - ta.nh.‘^{(3qimiZi))/3q2m2{l + Z 2 ) + 0{ml)] 


Zl,Z2 


= /3q2m2'^ Pp/2{zi\q){l - tanh^ {/3qimiZi) 
giving m 2 7^ 0 for 

Pm = 


(81) 


1 - E, PpnUh) tanhy/Sjimiz) 


(82) 


The theoretically predicted critical line (82) is in excellent agreement with the contour plot 
of m2 7^ 0 in the T — q 2 plane, shown in hg. 16 (left), computed numerically from (80). 
The plot shows that in the presence of 0{N"<) clones with higher numbers of receptors, the 
activation temperature of those with smaller promiscuity q 2 will deviate from the line T = q 2 . 


(iii) This is the case where deviations from the line T = q 2 are expected to be largest. From (41) 
we have 

mi = —( tanh /5(^^gimi-b N 

\ V ^ 

= ^ tanh {^{qimi + ^ qvC 


Pp- 2 {z\g^ tanh[/3(gimi(l + z) + ^2^20] 

Z 

Pp- 2 {z\ci^ tanh[/?gimi(l -b z)] + 0{N~^) 


(83) 
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Figure 16. Contour plot of m2 in the (T, (72) plane for q = (gi, ..., gi, (72, • ■ ■, 92) (left) and 
q = (91, ..., <7i, 92) (right). The blue dotted line represents the theoretical critical temperature 
line computed using respectively the self-consistent equations ( [80| ) (left) and (831 (right) together 
with the theoretical predictions for the critical temperature (|82[)(left) and (86((right). Deviations 
from the numerical results when (72 — 9i are due to the fact that (731 is obtained assuming 
7712 ^ rui, condition which is not satisfied when ^2 — 91. 


and 


7772 = ^ tanh 


(3{q2m2 + qvCmy) 





7F_i(2|q) taiih[/3(g2J«2 + qi'niiz)] 


(84) 


When clones with fewer receptors get active, we expect ttii = 0{1) and m 2 small. Taylor 
expanding for m 2 small and using the parity of Pn{z\q^ = Pn{—z\q) V n 

m 2 = Pp^i{z\q) [tanh(/jgir77iz) + (1 — ta.nh.'^{/3qimiz)(3q2m2) + 0{ml) 


- I3q2m2{l - ^2 7p-i(^|q) taiih^{l3qimiz)] + 0{m 


(85) 


we obtain for the activation temperature of the only clone with lowest number of receptors 

1 - Ez Pp-iiz\q) tanh.^{/3qimiz) 

in agreement with the contour plot of m 2 7^ 0 in the T 


( 86 ) 

q 2 found numerically, shown in fig. 
16 (right). Deviations are compatible with finite size effects and are largest for q 2 — qi where 


the assumption m 2 ttii is no longer valid. Monte Carlo simulations with 10 B clones, one 
of which has a lower promiscuity q 2 = 0.6, are shown in £g. IT at T = ^2 the pattern with 
lowest receptors is inactive and it is activated only at a much lower temperature T ~ 0.2. 
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Figure 17. Time evolution of B clones activation in Monte Carlo simulations with 10^ spins and 
q = (qi,..., < 71 , ( 72 ) (case (hi)) with qi = 0.8 (blue) and (72 = 0.6 (violet). The activation for the 
clone with less receptors (violet) decays to 0 for T = 0.5 (left) and stays non-zero for T = 0.2 
(right), i.e. at temperatures considerably lower than that in the absence of clonal interference 
T = q2. 


m 



Figure 18. Plot of m 2 as a function of T in the different cases: (1) q = (< 71 , ( 72 , • • • > Q 2 ), (2) 
q = (( 71 ,..., gi, ^ 2 , ■ • ■ 5 92 ), (3) q = (gi,..., gi, 52 )- Increasing the interference due to clones 
activated at a higher temperature, the m 2 intensity decreases. 


Finally in fig. 18 we plot m2 as a function of the temperature in the three different cases with 
promiscuity (i) q = (gi, ga, • • •, ^2), (h) q = (gi,..., gi, g2, • • •, g2 ), (hi) q = (gi,..., gi, q 2 )- 
The presence of clones with higher promiscuity gi and activation intensity mi, does not 
only affect the activation temperature of clones with lower promiscuity g = g2, but also the 
intensity m2 of their immune responses. In conclusion, clones with fewer triggered receptors 
are activated at a lower temperature and will produce a weaker response than clones with a 
higher number of active receptors. 
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5. Idiotypic interactions 


In this section we study the effect of the idiotypic interactions on the ability of the system to 
activate multiple clones in parallel. For simplicity, here we will assume homogeneous receptor 
promiscuities i.e. = c, V/i, so that the dynamical equations (18) read 


dm^ 

dt 






(87) 


pu 


Matrix A is positive definite for k G (0,1) and symmetric, hence the free-energy of the system is 
a Lyapunov function for the dynamics 123 and the system will converge to a steady state. Next, 
we compute the critical temperature at which clonal activation emerges in the steady state 

i\jj 

tanh(/jc (88) 

pu 

By summing over 

mx{Ar^)\pmf, = —(^ mxiA~^)xpC tanh(/3c CiA~^)p^mp))^ (89) 

fiX fiX pu 


using the inequality | tanh(a;)| < |a:| and averaging over the disorder we obtain 

'^mx{A~^)xpmp < (3c'^mx{A~'^)xpmp (90) 

piX pX 

Next we diagonalise the matrix A by means of the similarity transformation D = P^^AP, where 
D is the diagonal matrix constructed from the eigenvalues {p-p}^=i of A and P is the matrix of 
eigenvectors, which is unitary, i.e. P^^ = P^, since A is symmetric. Hence, we can rewrite the 
equation above for the transformed vector v = P~^m as 

- /3cv^D~V < 0 (91) 


which gives 

{/3c-fi.)>0 (92) 

u ^ ' 

yielding v = 0 for /3c < /Xmim where /i min = min^, [ly. The eigenvalues of A are 


= 1 - fc, deg{fii) = P/2 (93) 

P2 = 1 + deg{n2) = P/2 (94) 

hence, above the critical temperature Tc = c(/rmin)~^ = c/(l — k) we have m = 0, whereas below 

criticality m 7^ 0 is possible (and expected). Remarkably, the critical temperature increases with 

k, meaning that idiotypic interactions enhance immune system’s activation. Next, we investigate 
the structure of the states bifurcating away from m = 0 below Tf. and their stability. 
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Figure 19. Flow diagram in the plane {mi, m 2 ) (left) and Monte Carlo simulations with 10^ spins 
(right) for the dynamical system (95) with T = 1.7, c = 0.55, k = 0.7 (T^ = 1.83). 


5.1. Dynamical equations for two B clones 


As before, it is useful to consider first the toy model with P = 2 clones. For k ^ 0, the dynamics 
of different clones is always coupled, for all 7 > 0. We hrst look at the case 7 = 0 , where the 
system evolves according to the equations 


dmi 

dt 


dm 2 

dt 


(1 — c) tanh 


(3c 


1 - 


(mi + km2) ) + - 


tanh 


tanh 


(3c 


1 - k 


— mi 


(1 — c) tanh 


(mi - m 2 ) 

^^(fcmi + ms) j 


tanh 


Pc 

1-k 


Pc 

l-k 


(mi + m2) 


(mi + m2) 


— tanh 


Pc 

1 - k 


(mi - m2) 


- m 2 


(96) 


The state m = (0,0) is a fixed point of the dynamics at all temperatures, but we expect it to 
become unstable below = c/{l — k). In order to inspect the structure of the bifurcating state 
it is convenient to analyse the steady state equations in terms of the variables x = mi + m 2 and 
y = mi — m 2 
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Figure 20. Flow diagram in the plane {mi, m 2 ) (left) and Monte Carlo simulations with 10^ spins 
(right) for the dynamical system (95), for T = 0.65, c = 0.55,/c = 0.7 {Tc = 1.83). 


X = ctanh 
+ tanh 
y = ctanh 
— tanh 


(3cx 
1 - k 

(5c 


+ (l-c) 


tanh 


/3c 


1 - 

I3cy 


{x + {k — l)mi) 


l-k 

(5c 


+ (l-c) 


tanh 


l-F 


(5c 

l-F 


{x + {k - 1 ) 1712 ) 


{x + {k - 1 ) 1712 ) 


;{x + {k - l)mi) 


1 - P 

Assnming continnos bifurcations, we Taylor expand for small m near /5c = 1 — k i.e. at 
(5 = = 1 + e, obtaining, to leading orders 


X = c(5x — -c(5'^x'^ + (1 “ c)/3a; — 


1 — c 


-,(5^x 


(x + (fc - 1)7712)" 


3''' " ' 3(1 + A;)2 

+ (a; + (fc — l)77ii)^ — {x + {k — l)mi){x + (/c — 1 ) 7712 ) 


(96) 


y = c(5y- -cP^y^ + ^^ /37/ - {x + {k - 1 ) 1712 )^ 


1 

3' 


(l + k) 3(1 + A;)3 

+ {x + {k — l)77ii)^ — {x + {k — l)mi){x + {k — 1 ) 1172 ) ■ (97) 

We note that a; = 7/ = 0 is always a solution (corresponding to m = (0, 0)). A solution 7/ 7^ 0 is 
not possible as in (97) terms O(e^) do not simplify. Hence, 7/ = 0 is the only solution, implying 
777i = 7772 . Ih coutrast, in (96) first order terms simplify, hence we can have a: 7^ 0. Therefore, 
mixtures bifurcate from m = (0,0) in a symmetric fashion m = 777(1,1). We can compute the 



































amplitude m, Taylor expanding (95) at the steady state near Tc, for small m, 


m = /3m — 




yielding, at /3 = 1 + e, 


m? = 


3e 


( 98 ) 


(99) 


l + 3c 

or the trivial solntion m = 0. Flow diagrams and Monte Carlo simnlations are in agreement with 
theoretical predictions, showing that close to criticality both B clones are activated with the same 
intensity (£g. 19). Next we investigate the stability of symmetric solntions by compnting the 
eigenvalnes of the Jacobian 


j {IV 


dF, 


(1), 


m 


drrii, 


= — (^^tanh(/3c^^^(A ^)pxmp))^ - trip ,(100) 

p\ 


One has 

Jp, = l3c{A-^)p,{l - (tanh2(/3c5^e\A-i),Amp)^) - 5p, + 0{N-^) 

pX 

and snbstitnting m* = m(l, 1) we hnd for the eigenvalnes 


Ai = 


Ao — 




l-k 
/3c(l — c) 
1 + k 


l-k 
1 — tanh 


m(3c 
1 - k 


(S)-» 

/3c2 


+ 


1 - k 


l-k 
- 1. 


tanh 


/ 2m/3c 
\l-k 


These can be calcnlated analytically near criticality i.e. at /3c/l — k = 1 + e where 
—2e — 6ec 

Ai = 

A2 = 


1 + 3c ’ 

—2e + 6ec — 2k + 2ke — Akc — Akce 


(1 + 3c) (1 + k) 


and for T —)■ 0, where 

Ai = —1 , 
I3c^ 


Ao — 


l-k 


- 1 


( 101 ) 

( 102 ) 

(103) 

(104) 

(105) 

(106) 
(107) 


We dednce that Ai is always negative, hence the stability of symmetric solntions is determined by 
A2. A plot of A2 as a fnnction of T is shown in £g. 21 (left) for hxed dilntion and B-B interaction 
strength k. The critical line where A2 = 0 in the T — c plane for different values of k is shown 
in £g. 21 (right). Remarkably, the region (S) where activation of parallel immune responses is 
accomplished with the same intensity gets wider as k increases. For 7 > 0, the system evolves 
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Figure 21. Left: A 2 T as a function of T for c = 0.3 and k=0.2. The value at T = 0 matches the 
analytical prediction A 2 T = c^/1 — k (red marker). Right: Contour plot of A 2 = 0 with different 
interaction strengths k = 0.8, k = 0.5, fc = 0.2 in the parameter space {T,c). The region where 
clones are activated in parallel with the same intensity (S) becomes wider as k increases. 


according to the equations 


dmi 

dt 


tanh 


dm2 

dt 


tanh 


l3c 


1 - 

jdc 




(mi + km2) 
{kirti + m 2 ) 


— mi 


- m2 


(108) 


One can show again that m = (0, 0) is the only hxed point above Tc = c/[1 — k), and symmetric 
mixtures m = m(l, 1) bifurcate away from m = (0,0) at Tc- Now, however, symmetric mixtures 
remain stable for any T < (see Sec. 5.2.1), with the intensity m of the symmetric activation 
found from 


“ = tanh 


(109) 


Hence for 7 > 0 the system reduces to independent Curie-Weiss ferromagnets, even in the presence 
of idiotypic interactions, with the critical temperature Tc = c/(l — k) increasing with the strength 
of interactions. 


5.2. Generalisation to P clones 

In this section we turn to the case of P = , 0 < 5 < 1 B clones. Again, we expect clones to 

be activated at Tc = c/(l — k) in a symmetric fashion. For 5 < 7 we expect symmetric solutions 

to be stable for all T < Tc, whereas for 5 > 7 we expect them to destabilise at low temperature. 

Without loss of generality, we can set P = aN'^ with a = 1 for 5 = 7 and the cases 5 < 7 and 

(5 > 7 retrieved in the limits a —)■ 0 and a —)■ 00 respectively. 

To inspect the behavior of the system near criticality, it is convenient to write the steady 
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state equations in terms of the rescaled matrix A = A/{1 — k), with eigenvalues 

degini) = P/2 

deg{ii2) = P /2 . 


hi = 1 , 

1 + k 

h 2 = 


1 - k 


( 110 ) 

( 111 ) 


Taylor expanding the steady state equations (88) for small m near Tc, we have, with /3 = 
N'T 


pu 


/33 ATT' 

ir~ 




Pi,’^1,P2,1^2,P3,1^3 


( 112 ) 


and averaging over the disorder gives 
= ^ J2^A)-^rnp - y 

p p 

To make progress, we project on the eigenvectors of A“^; 


"^p - J2^A)^^mp 

p T7T 


pp "'p 




( 113 ) 


J = 1 




ajV'fi' = 


i=l i,j,k 


( 114 ) 

P7^P j,k 

Next we split the sum into j = 1,..., P/2 and j = P/2 + 1,..., P, substitute the eigenvalues 
(110), (111), and equating linear terms we get 

-1^ p . , j, 


^o_i!i“(l-/3)+ ^ a,!iO> 1 - = 0 


i=i 


j=P/2+l 


1+v “ 

At /3 = 1 linear terms cancel only if Yl,^=p/ 2 +i = 0 V /i, which implies vj/^ = 0 V j = 

^ + 1,..., P, V p. We note that the eigenvectors of A~^ are in the form 


= ^i^p + mp+p/ 2 ) , j = 1, ..,P/2 , 

= ^{mp- mp+Pi2) , j = P/2 + 1,.., P , 

,(i) 


(116) 

(117) 


Therefore we have, for all p, rrip = m^+p /2 and Vp = rrip V j = l,...,P/2. This implies 
nip = (^j'kiT'pi hence = 1- these conditions into (114) 

/33 


mp{-l +--ml-= ^ , 

p¥=p 


we have at /3 = 1 + e. 


nine -mf, — ——mu|mk = 0 , 

P 2 ^ N'T ' 


( 118 ) 


( 119 ) 


31 










where |mp = Ylp'^% yielding = 0 


or 


3c 


2 o I |2 

m,, = 3e-m . 

/i I I 


( 120 ) 


Each clone activation rrip depends on the whole vector m, hence all clones will have the same 
activation strength rrip, G (—m, +m, 0). Assuming the non-zero components are a fraction (j) < a 
of the total number of components P, we have |mp = yielding 

^ - ( 121 ) 

1 + 3ca0 

We will see in the next section that while for fc = 0 stability of symmetric mixtures is only ensured 
for 0 = 1, in the presence of idiotypic interactions, he. for fc 7^ 0, values (j) ^ 1 are possible. 


5.2.1. Linear stability analysis and phase diagram In this section we study the stability of 
symmetric clonal activation, m = m(l,..., 1 , 0 ,... 0) with n = acpN"' activated clones, below 
criticality, in the presence of idiotypic interactions. To this purpose, we study the eigenvalues of 


the Jacobian of the dynamical system (87), which has a block structure with diagonal elements 
given for /i < n by 

Jpp = /3(A“^)^^(1 - (tanh^[/?(^(A~^)p^mp + ^ ^^(A~^)pAmp)])^) - 1, (122) 

P P,A7:^ 

and for /i > n by 

Jpp = /5(A~^)/,;,[1 - (tanh^/3( ^ ^^(A^^)pAmp)])^ - 1, (123) 


p>Vm 


while off-diagonal elements are, for fi,i> < n 


JflU 




(tanh^[/3( ^^(A ^)pxmp))^ - 1] 

p,\^p 


(124) 


and null otherwise. For A —)■ 00 the matrix becomes diagonal and eigenvalues are given by the 
diagonal elements. At the symmetric hxed point, using X]p(A^^)pA = 1/(1 — A;) V A, we get 


Ai = 


Ao — 




l + k 

l + k 


1 — Pn{z) tanh^(/3(m(l + 2)) j — 1 deg{Xi) = n , 


(125) 


1 — Y^ Pn{z) ta.nh^0mz) 1 — 1 


deg{\i) = P — n (126) 


with Pn{z) dehned by (77) with q = c(l,..., 1) and given by Pn{z) = e where Iz{x) 


is the modihed Bessel function of the hrst kind. Close to criticality, by Taylor expanding in power 
of e = /j — 1 we obtain 

—2e — A:(l + Scacp) 


Ai = 
A2 = 


(1 -|- A;)(l -|- Scacj)) 
e — k{l + 3ca(j)) 
(1 J- A:)(l J- Scacj)) 


(127) 

(128) 
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In contrast to the case fc = 0 where A2 > 0 and symmetric activation is stable only if it involves 
all clones, in the presence of idiotypic interactions, i.e. for k ^ 0, both eigenvalues are negative at 
criticality, showing the emergence of local minima where not all the clones are activated. In the 
opposite limit, i.e. T —0 we get 


Ai = 


Ao — 




1 + k 
1 + k 




(129) 


(130) 


For (5 < 7, a = 0, hence using the properties of Bessel functions /i(0) = 0 and /o(0) = 1, we get 
Ai < 0 and A2 > 0. Since A2 has degeneracy P — n, symmetric mixtures m = m(l,..., 1) i.e. with 
n = P will be stable for all T < T^.. In contrast, for 5 = 7 he. a 7^ 0, Ai > 0 as Ii{x) > 0 for any 
X > 0, meaning that symmetric mixtures are unstable at low temperature for any n. Finally, for 
(5 > 7 one has a —)■ cxd, in the thermodynamic limit so In{o:(j)c) ~ (q!(/)c)“^/^ V u, and symmetric 
mixtures will gain stability at low temperature as iV —)■ cx). 

Let us now compute the critical line in the phase diagram where symmetric mixtures become 
unstable. We note that A2 > Ai for T —)■ 0 and for T —)■ Tc, and deduce that A2 > Ai for all 
T < Tc. Hence, stability of n-mixtures m = m(l,..., 1, 0,..., 0) is given by the region where 
A2 < 0. In £g. 22 (left) we show the critical lines where A2 gets zero in the space of scaled 


parameters T = T/c, (f) = (j)c for different values of k and a = 1. As k increases the region where 
clones are activated with the same intensity gets wider. When A2 destabilises, symmetric mixtures 
can only be stable for n = P, i.e. for 0 = 1, in the region where Ai < 0. A contour plot of Ai = 0 
for 0 = 1 in the T — k plane is shown in fig. 22 (right). 



- k=0.05 


- k=0.5 



Figure 22. Left: Contour plot of A 2 = 0 (obtained from (126l) for a = 1, as a function of the 
scaled parameters [T = T/c^(j) = (jic, for different value of B-B interaction strength k. Increasing 
k the region where symmetric mixtures are stable (to the right of the critical line) becomes wider. 
Right: Contour plot of Ai = 0 (obtained from ( 125| )) for a = 0 = 1, as a function of the scaled 
temperature T = T/c and strength k of the idiotypic interactions. 
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6. Antigen effect 


In this section we investigate the effect of antigens on the basal activity of the immune system 
analysed in the previous sections. We will carry out the analysis for homogeneous promiscuities 
and in the absence of idiotypic interactions. We suppose to have A = antigens, in the presence 
of which dynamical equations ([T^ become 

dru ^r 

= —(^^tanh[/3^^''(m^ + , (131) 

v=l 


where j3 = jSc and ^lJy is the antigenic held associated to B clone z/. At the steady state we have 
two sets of equations, those for clones /i = 1 ,A complementary to the incoming viruses and 
those for non-activated clones u = A + 1,..., P, performing basal activities: 

A p 

rrif, = (tanh(/3(m^ + + y^^rn^C))^ , /i = l,...,A, (132) 

A P 

rriu = {ta.nh0 {mu + , p = A + l,...,P. (133) 

pj^p p^p 

We hrst look at the case A <C iV'’', where the equations can be written as 


mu 


(izP{z\m) tanh(/3(m^ + -0^ + z)) , /i = 1,..., A , 
<lzP{z\m) tanh(/3(m,y + '^)) ) z/ = A + 1,..., P , 


(134) 

(135) 


where P{z\Yn) is the large N limit of P^(2;|{mp, 1}) which was dehned in (44), and only depends 
on the vector m of basal activation. 

For P -C , we have P(2;|m) = {5{z — J2u — '^(^) equations decouple, 

hence there is no interference between infected and non-infected clones 


rup = tanh(/3(mp +-z/zp)) , /i = l,...,A, (136) 

mu = tanh(/3mi.) , z/ = A + l,...,P. (137) 

For the infected clones, the presence of the held induces hysteresis ehect in the clonal activation 
[23] . This may explain immunological memory ehects |2il [12], without the requirement of 
dedicated memory cells: after an infection, the responsive B cells, may retain a non-zero activation 
as the antigen is fought and its concentration is sent to zero, and on a successive encounter with 
the same antigen they will provide a higher and faster response. 

For P = 0{N'^), P(2:|m) in (134), (135), has a hnite width, and both antigen-induced and 
basal activities reduce, due to clonal interference, and hysteresis cycles become smaller. However, 
the presence of antigens does not ahect the basal activity of uninfected clones u = A + 1,..., P, 
as long as A N'^. 

Next we investigate the scenario where the number of antigens is A = 0{N"^) and ask whether 
the system is able to hght against all of them in parallel. For simplicity, we set P = N'^ and 
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A = (piN'^, with 01 denoting the fraction of infected clones and we assnme that all viruses have 
the same concentrations i.e. xj) = 0(1,..., 1, 0,..., 0). This leads to the steady state equations 

A p 

rUf, = {tanh0{mf, + xlj+ + i>)C + ^ = 

u=l u=A+l 

A P 

rUj, = (tanh(/3(mj, + ^(rUp + 0)^^ + ^ u = A + 1,... ,P 

p=l p=A+l 

(138) 

Now antigen interference on the basal activity is relevant and will affect the activation of the 
non-infected clones. In the small field limit, we can Taylor expand (138) near 0 = 1 + e and small 
m^, obtaining 


1 / 

Trip ~ (1+ e)m^ + 0^(l+ e) - - ((m^ + 0^)^ + ^(m^ + 0^) ^(rup + 0p) 


For fx < A (infected clones) we have for 0 <C e 

—0(1 + e) T A 

mp = - =mi, ia = 1,...,A. 


(139) 


(140) 


where the expansion holds for held 0 e, otherwise rrip = 0{1). For fi > A (0p = 0) we have 

mp ~ (1 + e)mp - ^ (ml + ^^(mp + xPpfj (141) 


hence = 0 is always a solution (uninfected clones may not be activated) together with 


3c 


^ J2^rnp + 0p)2 = 


p 


= 3e-— ( + + > n = A + 1,...,P . 

p=A+l p=l 

This shows that uninfected clones are symmetrically activated, each one with intensity 


m2 = 


3e - 3 c.)'iVj^( 2(1 + 6)/e + 1) 3^(1 - 2c4>X^Y) 


(142) 


(143) 


1 + 3c02 1 T 3c02 

where 02 < 1 ~ 0i is the fraction of active uninfected clones. Hence, for 0 <C e, clonal activation 
close to criticality will have the form m* = (mi,..., mi, m 2 ,..., m 2 , 0,..., 0). However, upon 
increasing the fraction 0i of infected clones or the antigenic held 0, (143) shows that non¬ 
zero values of m2 may become impossible and uninfected clones may get activated at a lower 
temperature (similarly to clones with smaller numbers of triggered receptors we dealt with in Sec. 
1^. This is conhrmed by numerical results in Fig. 23, showing that the response of uninfected B 
clones and their activation temperature decrease for increasing fractions of infected clones. This 
results in a reduced basal activity of the immune system, which is vital to keep cells signaled 
and accomplish homeostasis [IH]. In hg. 24 we study the impact of antigen concentration on 
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Figure 23. 3D plot of the activation m 2 of non-infected clones versus the fraction of infected 
clones 4>i = (/)ic and the scaled temperature T = T/c^ for fixed ip = 0.1. Increasing <f>i both the 
intensity and the critical activation temperature decrease, due to the antigenic interference. 



-m - l/r=0.01 


(if=0.05 


ifr=Q.2 


-A- ifr=Q.5 


Figure 24. Critical line for the activation of uninfected clones in the space of scaled parameters 
T = T/c, (j) = (pc, otained from the condition m 2 ^ 0 (141), for different values of the antigenic 
field Ip. Increasing ip, the region where uninfected clones are signaled shrinks. 


the critical temperature at which uninfected B clones become responsive, by plotting the critical 
temperature versus the fraction of infected clones, for different values of antigen concentration. 
This shows that as the fraction of infected clones and the held increase, the basal activity is more 
and more compromised. 

Next, we inspect the stability region of m* = (mi,..., mi, m 2 ,..., m 2 , 0,..., 0) by looking at 
the eigenvalues of the Jacobian of the dynamical system (131). This has a diagonal structure, in 
the thermodynamic limit, where off-diagonal elements become negligible and diagonal terms are 
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n2 


u=ni-\-l 


ni 


u=l 


for n <ni = 

=/3(1 - (tanh^(^(m^ +-0 + ^ + 'ip)C')))fi) 

for 77-1 < /r < n 2 = 02 -^"’' 

=/5(1 - (tanh^(/3(m^ + ^ + V’)D))^) 

and for fi > n 2 

J'/.M = /3(1 - (tanh^(^( ^ + ^)D))^) 




u=ni+l 


n\ 


u=l 


n2 


n\ 


(144) 


(145) 


(146) 


v=ni+l 


u=l 


Evaluating the Jacobian at the symmetric hxed point m* and introducing the distribution 
Pn{z) = {6{z - Y!^^i D) , V £, gives 

Ai = /3{l-y^Pni(^:i)-Pn2(^2) tanh^[/3((7/> + mi)(l + zi) + m 2 ^ 2 )]}-l deg{Xi)=ni (147) 

^ 1 , 2:2 

X2 = ^{l-y^Pniizi)Pn2iz2)tanh^[/3{{i> + mi)zi + 7772(1 + 2 : 2 )]} -1 deg{X2)=n2 (148) 

^ 1 , 2:2 

A 3 = /3{l-y^Pni(2:i)Pn2(^2)tanh^[^((77 + r?7i)^i + 77722 : 2 )]} -1 deg{X‘i)=P-ni-n 2 (149) 


•^1,22 


with mi,m 2 following from ( |140 ), (141) as 

7771 = E ^ni(+)^412(2:2) tanh0{{^jJ + 777i)(l + 2:1) + 7772^2)] 


- 21,-22 


m 2 ^ ^ Pni(zi)Pn 2 (z 2 ) taah[l3({il2 + m 2 )Z] + 7712(1 + 1 : 2 ))] 


(150) 


- 21,-22 


and Pn{z) = e^‘^Iz{4>c) for n = (pN"', where Iz{x) is the modified Bessel function of the first kind. 
In Fig. 25 (left) we show the critical lines Ai = 0, A 2 = 0 and A 3 = 0 in the space of scaled 
parameters T = T/c, 02 = (f> 2 C. As temperature is lowered, A 3 is the hrst eigenvalue to destabilise, 
meaning that clonal activation will get in the form m = ( 7771 ,... , 7771 , 7772 ,... , 7772 ). Decreasing 
the temperature further, the system will first prioritise activation of infected clones, while keeping 
activation of uninfected clones symmetric, and later, at low temperature, will activate uninfected 
clones in a hierarchical fashion (see fig. 25, right panel). Near the critical temperature symmetric 
mixtures with 02 < 1 — 0i are stable. To investigate the optimal value of 02 we calculate the 
free-energy as a function of 02 for hxed T, 0i 

P( 01 , 02 ) = “W X] Pni( 2 :i)Pn 2 ( 2 : 2 ) log (2 COsh( 0 c( 777 i + 'ljj)zi + 77722 : 2 )) + -( 0 l 777 i + (p 2 ml) (151) 

21 , 2:2 

which is minimal (hg. [26| at 02 = 1 — 0 i, meaning that the system will keep all uninfected clones 
signalled. 
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Figure 25. Left: Phase diagram in the space of scaled parameters f = T/c, (^2 = </> 2 C with 
(jji = 0.4. Lines represent contours of Ai = 0 (circles), A 2 = 0 (triangles) and A 3 = 0 (squares). 
To the right of the line A 3 = 0 solutions where uninfected clones are partially activated are 
stable. Lowering down the temperature and crossing the line Ai = 0 the mi symmetric mixtures 
destabilise, meaning that infected clones are hierarchically activated. Crossing the line A 2 = 0 
the 7712 symmetric mixtures destabilise, and uninfected clones are hierarchically activated. Right: 
Plots of Ai = 0 (circles), A 2 = 0 (triangles) in the space T = T/c, ^2 = 4>2C for ^2 = 1 — <(>i- 


t(02) 



Figure 26. Free energy as a function of the fraction of active uninfected clones (^ 2 ) for (pi = 0.5, 
T = 0.7,c = 0.2,'i/; = 0.05. 


7. Conclusions 

In this paper we presented a minimal model for interacting cells in the adaptive immune systems, 
constituting of B cells, T cells and antigens. Our model is able to capture important collective 
features of the real immune system, such as the ability of simultaneously handling multiple 
infections, the dependence of B clones’ activation on the number of receptors on clonal surface and 
the role of idiotypic interactions in enhancing parallel response to multiple infections. We analysed 
the dynamics of the system’s order parameters quantifying the B clones’ activation via linear 
stability analysis and Monte Carlo simulations. We found the regions, in the parameters’ space, 
where the system activate B clones symmetrically, i.e. with the same intensity, or hierarchically, 
whereby the system prioritises responses to particular infections. We showed that clones with 
fewer receptors are less likely to be activated (Sec. and idiotypic interactions contribute to 
the overall stability of the immune system, preventing unwanted activation and increasing the 
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region where all clones are equally activated and ready to start an immune response upon arrival 
of new infections (Sec. |^. Furthermore, we investigated how the immune system responds to 
antigens, showing in particular that multiple antigens create an interference that leads to less 
effective response to individual antigens and a reduction in the basal activity of the uninfected B 
cells (Sec. [^. For higher noise level the system tends to simultaneously hght all the antigens, 
despite losing in terms of response strength, whereas for lower noise level it will prioritise some 
infections over others. Finally, we showed that short term memory may emerge as an hysteresis 
effect, without the requirement of dedicated memory cells. 
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